{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d068cb54",
   "metadata": {},
   "source": [
    "# Statistics and inference for one and two sample Poisson rates\n",
    "\n",
    "Author: Josef Perktold\n",
    "\n",
    "This notebook provides a brief overview of hypothesis tests, confidence intervals and other statistics for Poisson rates in one and two sample case. See docstrings for more options and additional details.\n",
    "\n",
    "All functions in `statsmodels.stats.rates` take summary statistics of the data as arguments. Those are counts of events and number of observations or total exposure. Some functions for Poisson have an option for excess dispersion. Functions for negative binomial, NB2, require the dispersion parameter. Excess dispersion and dispersion parameter need to be provided by the user and can be estimated from the original data with GLM-Poisson and discrete NegativeBinomial model, respectively.\n",
    "\n",
    "Note, some parts are still experimental and will likely change, some features are still missing and will be added in future versions.\n",
    "\n",
    "[One sample functions](#One-sample-functions)  \n",
    "[Two sample functions](#Two-sample-functions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59f2450e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.stats.rates as smr\n",
    "from numpy.testing import assert_allclose\n",
    "from statsmodels.stats.rates import (  # functions for 1 sample; functions for 2 sample; power functions; list of statistical methods\n",
    "    confint_poisson,\n",
    "    confint_poisson_2indep,\n",
    "    confint_quantile_poisson,\n",
    "    etest_poisson_2indep,\n",
    "    method_names_poisson_1samp,\n",
    "    method_names_poisson_2indep,\n",
    "    nonequivalence_poisson_2indep,\n",
    "    power_equivalence_neginb_2indep,\n",
    "    power_equivalence_poisson_2indep,\n",
    "    power_negbin_ratio_2indep,\n",
    "    power_poisson_diff_2indep,\n",
    "    power_poisson_ratio_2indep,\n",
    "    test_poisson,\n",
    "    test_poisson_2indep,\n",
    "    tolerance_int_poisson,\n",
    "    tost_poisson_2indep,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cab009f6",
   "metadata": {},
   "source": [
    "## One sample functions\n",
    "\n",
    "The main functions for one sample Poisson rates currently are test_poisson and confint_poisson. Both have several methods available, most of them are consistent between hypothesis test and confidence interval. Two additional functions are available for tolerance intervals and for confidence intervals of quantiles.  \n",
    "See docstrings for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "199276f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "count1, n1 = 60, 514.775\n",
    "count1 / n1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c57f1dfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_poisson(count1, n1, value=0.1, method=\"midp-c\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae000bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "confint_poisson(count1, n1, method=\"midp-c\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d27cff72",
   "metadata": {},
   "source": [
    "The available methods for hypothesis tests and confidence interval are available in the dictionary `method_names_poisson_1samp`.  See docstring for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5b89f88",
   "metadata": {},
   "outputs": [],
   "source": [
    "method_names_poisson_1samp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac9d09a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "for meth in method_names_poisson_1samp[\"test\"]:\n",
    "    tst = test_poisson(count1, n1, method=meth, value=0.1, alternative=\"two-sided\")\n",
    "    print(\"%-12s\" % meth, tst.pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f224485",
   "metadata": {},
   "outputs": [],
   "source": [
    "for meth in method_names_poisson_1samp[\"confint\"]:\n",
    "    tst = confint_poisson(count1, n1, method=meth)\n",
    "    print(\"%-12s\" % meth, tst)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31a7f349",
   "metadata": {},
   "source": [
    "Two additional functions are currently available for one sample poisson rates, `tolerance_int_poisson` for tolerance intervals and `confint_quantile_poisson` for confidence intervals of Poisson quantiles. \n",
    "\n",
    "Tolerance intervals are similar to prediction intervals that combine the randomness of a new observation and uncertainty about the estimated Poisson rate. If the rate were known, then we can compute a Poisson interval for a new observation using the inverse cdf at the given rate. The tolerance interval adds uncertainty about the rate by using the confidence interval for the rate estimate.\n",
    "\n",
    "A tolerance interval is specified by two probabilities, `prob` is the coverage of the Poisson interval, `alpha` is the confidence level for the confidence interval of the rate estimate.  \n",
    "Note, that probabilities cannot be exactly equal to the nominal probabilites because counts are discrete random variables. The properties of the intervals are specified in term of inequalities, coverage is at least `prob`, coverage of the confidence interval of the estimated rate is at least `1 - alpha`. However, most methods will not guarantee that the coverage inequalities hold in small samples even if the distribution is correctly specified.\n",
    "\n",
    "In the following example, we can expect to observe between 4 and 23 events if the total exposure or number of observations is 100, at given coverage `prob` and confidence level `alpha`. The tolerance interval is larger than the Poisson interval at the observed rate, (5, 19), because the tolerance interval takes uncertainty about the parameter estimate into account."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05f9882e",
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure_new = 100\n",
    "tolerance_int_poisson(\n",
    "    count1,\n",
    "    n1,\n",
    "    prob=0.95,\n",
    "    exposure_new=exposure_new,\n",
    "    method=\"score\",\n",
    "    alpha=0.05,\n",
    "    alternative=\"two-sided\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97342d26",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "stats.poisson.interval(0.95, count1 / n1 * exposure_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86de470d",
   "metadata": {},
   "source": [
    "Aside: We can force the tolerance interval to ignore parameter uncertainty by specifying `alpha=1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45696c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance_int_poisson(\n",
    "    count1,\n",
    "    n1,\n",
    "    prob=0.95,\n",
    "    exposure_new=exposure_new,\n",
    "    method=\"score\",\n",
    "    alpha=1,\n",
    "    alternative=\"two-sided\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6e62fdf",
   "metadata": {},
   "source": [
    "The last function returns a confidence interval for a Poisson quantile. A quantile is the inverse of the cdf function, named `ppf` in scipy.stats distributions.\n",
    "\n",
    "The following example shows the confidence interval for the upper bound of the Poisson interval at cdf probability 0.975. The upper confidence limit using the one-tail coverage probability is the same as the upper limit of the tolerance interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a7df495",
   "metadata": {},
   "outputs": [],
   "source": [
    "confint_quantile_poisson(\n",
    "    count1,\n",
    "    n1,\n",
    "    prob=0.975,\n",
    "    exposure_new=100,\n",
    "    method=\"score\",\n",
    "    alpha=0.05,\n",
    "    alternative=\"two-sided\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5107acc",
   "metadata": {},
   "source": [
    "## Two sample functions\n",
    "\n",
    "Statistical function for two samples can compare the rates by either the ratio or the difference. Default is comparing the rates ratio.\n",
    "\n",
    "The `etest` functions can be directly accessed through `test_poisson_2indep`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6c7d341",
   "metadata": {},
   "outputs": [],
   "source": [
    "count1, n1, count2, n2 = 60, 514.775, 30, 543.087"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aba830a",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_poisson_2indep(count1, n1, count2, n2, method=\"etest-score\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52d1cd76",
   "metadata": {},
   "outputs": [],
   "source": [
    "confint_poisson_2indep(count1, n1, count2, n2, method=\"score\", compare=\"ratio\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3642594e",
   "metadata": {},
   "outputs": [],
   "source": [
    "confint_poisson_2indep(count1, n1, count2, n2, method=\"score\", compare=\"diff\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b30dfd37",
   "metadata": {},
   "source": [
    "The two sample test function, `test_poisson_2indep`, has a `value` option to specify null hypothesis that do not specify equality. This is useful for superiority and noninferiority testing with one-sided alternatives.\n",
    "\n",
    "As an example, the following test tests the two-sided null hypothesis that the rates ratio is 2. The pvalue for this hypothesis is 0.81 and we cannot reject that the first rate is twice the second rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85046eb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_poisson_2indep(count1, n1, count2, n2, value=2, method=\"etest-score\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4831cf3f",
   "metadata": {},
   "source": [
    "The `method_names_poisson_2indep` dictionary shows which methods are available when comparing two samples by either rates ratio or rates difference.\n",
    "\n",
    "We can use the dictionary to compute p-values and confidence intervals using all available methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6a6aec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "method_names_poisson_2indep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "302b0bf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "for compare in [\"ratio\", \"diff\"]:\n",
    "    print(compare)\n",
    "    for meth in method_names_poisson_2indep[\"test\"][compare]:\n",
    "        tst = test_poisson_2indep(\n",
    "            count1,\n",
    "            n1,\n",
    "            count2,\n",
    "            n2,\n",
    "            value=None,\n",
    "            method=meth,\n",
    "            compare=compare,\n",
    "            alternative=\"two-sided\",\n",
    "        )\n",
    "        print(\"   %-12s\" % meth, tst.pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b306fa76",
   "metadata": {},
   "source": [
    "In a similar way we can compute confidence intervals for the rate ratio and rate difference for all currently available methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2658911",
   "metadata": {},
   "outputs": [],
   "source": [
    "for compare in [\"ratio\", \"diff\"]:\n",
    "    print(compare)\n",
    "    for meth in method_names_poisson_2indep[\"confint\"][compare]:\n",
    "        ci = confint_poisson_2indep(\n",
    "            count1, n1, count2, n2, method=meth, compare=compare\n",
    "        )\n",
    "        print(\"   %-12s\" % meth, ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65ce8685",
   "metadata": {},
   "source": [
    "We have two additional functions for hypothesis tests that specify interval hypothesis, \n",
    "`tost_poisson_2indep` and `nonequivalence_poisson_2indep`.\n",
    "\n",
    "The `TOST` function implements equivalence tests where the alternative hypothesis specifies that the two rates are within an interval of each other.\n",
    "\n",
    "The `nonequivalence` tests implements a test where the alternative specifies that the two rates differ by at least a given nonzero value. This is also often called a minimum effect test. This test uses two one-sided tests similar to TOST however with null and alternative hypothesis reversed compared to the equivalence test.\n",
    "\n",
    "Both functions delegate to `test_poisson_2indep` and, therefore, the same `method` options are available."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74e1de07",
   "metadata": {},
   "source": [
    "The following equivalence test specifies the alternative hypothesis that the rate ratio is between 0.8 and 1/0.8. The observed rate ratio is 0.89. The pvalue is 0.107 and we cannot reject the null hypothesis in favor of the alternative hypothesis that the two rates are equivalent at the given margins. Thus the hypothesis test does not provide evidence that the two rates are equivalent.\n",
    "\n",
    "In the second example we test equivalence in the rate difference, where equivalence is defined by margins (-0.04, 0.04). The pvalue is around 0.2 and the test does not support that the two rates are equivalent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ee7a0ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "low = 0.8\n",
    "upp = 1 / low\n",
    "\n",
    "count1, n1, count2, n2 = 200, 1000, 450, 2000\n",
    "\n",
    "tost_poisson_2indep(count1, n1, count2, n2, low, upp, method=\"score\", compare=\"ratio\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d92005e",
   "metadata": {},
   "outputs": [],
   "source": [
    "upp = 0.04\n",
    "low = -upp\n",
    "tost_poisson_2indep(count1, n1, count2, n2, low, upp, method=\"score\", compare=\"diff\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7329c90",
   "metadata": {},
   "source": [
    "The function `nonequivalence_poisson_2indep` tests the alternative hypothesis that the two rates differ by a non-neglibile amount.\n",
    "\n",
    "In the following example, the alternative hypothesis specifies that the rate ratio is outside the interval (0.95, 1/0.95). The null hypothesis is that the ratio ratio is in the interval. If the test rejects the null hypothesis, then it provides evidence that the rate ratio differ by more than the unimportant amount specified by the interval limits.\n",
    "\n",
    "A note on the relationship between point hypothesis test and interval hypothesis test in large samples. The point null hypothesis of test_poisson_2indep will reject any small deviation from the null hypothesis if the null hypothesis does not hold exactly and the sample size is large enough. The nonequivalence or minimum effect test will not reject the null hypothesis in large samples (sample approaches infinite) if rates differ by not more than the specified neglibible amount.\n",
    "\n",
    "In the example neither the point nor the interval null hypothesis are rejected. We do not have enough evidence to say that the rates are statistically different.\n",
    "Following that, we increase the sample size 20 times while keeping observed rates constant. In this case, the point null hypothesis test is rejected, the pvalue is 0.01, while the interval null hypothesis is not rejected, the pvalue is equal to 1.\n",
    "\n",
    "Note: The nonequivalence test is in general conservative, its size is bounded by alpha, but in the large sample limit with fixed nonequivalence margins the size approaches alpha / 2. If the nonequivalence interval shrinks to a single point, then the nonequivalence test is the same as the point hypothesis test. (see docstring)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32ad0255",
   "metadata": {},
   "outputs": [],
   "source": [
    "count1, n1, count2, n2 = 200, 1000, 420, 2000\n",
    "low = 0.95\n",
    "upp = 1 / low\n",
    "nf = 1\n",
    "nonequivalence_poisson_2indep(\n",
    "    count1 * nf,\n",
    "    n1 * nf,\n",
    "    count2 * nf,\n",
    "    n2 * nf,\n",
    "    low,\n",
    "    upp,\n",
    "    method=\"score\",\n",
    "    compare=\"ratio\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61beae8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_poisson_2indep(\n",
    "    count1 * nf, n1 * nf, count2 * nf, n2 * nf, method=\"score\", compare=\"ratio\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62b98f47",
   "metadata": {},
   "outputs": [],
   "source": [
    "nf = 20\n",
    "nonequivalence_poisson_2indep(\n",
    "    count1 * nf,\n",
    "    n1 * nf,\n",
    "    count2 * nf,\n",
    "    n2 * nf,\n",
    "    low,\n",
    "    upp,\n",
    "    method=\"score\",\n",
    "    compare=\"ratio\",\n",
    ").pvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "212ca284",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_poisson_2indep(\n",
    "    count1 * nf, n1 * nf, count2 * nf, n2 * nf, method=\"score\", compare=\"ratio\"\n",
    ").pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dc3f3bc",
   "metadata": {},
   "source": [
    "## Power\n",
    "\n",
    "Statsmodels has limited support for computing statistical power for the comparison of 2 sample Poisson and Negative Binomial rates. Those are based on Zhu and Lakkis and Zhu for ratio comparisons for both distributions, and basic normal based comparison for the Poisson rate difference. Other methods that correspond more closely to the available methods in the hypothesis test function, especially Gu, are not yet available.\n",
    "\n",
    "The available functions are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2479efb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "power_poisson_ratio_2indep\n",
    "power_equivalence_poisson_2indep\n",
    "power_negbin_ratio_2indep\n",
    "power_equivalence_neginb_2indep\n",
    "\n",
    "power_poisson_diff_2indep"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
